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Abstract 

The focus of this research is on the development of analysis and sensitivity 
analysis equations for nonlinear, transient heat transfer problems modeled by p- 
version, time discontinuous finite element approximation. The resulting matrix 
equation of the state equation is simply in the form oiA(x)x = c , representing a single 
step, time marching scheme. The Newton-Raphson’s method is used to solve the 
nonlinear equation. Examples are first provided to demonstrate the accuracy 
characteristics of the resultant finite element approximation. A direct differentiation 
approach is then used to compute the thermal sensitivities of a nonlinear heat transfer 
problem. The report shows that only minimal coding effort is required to enhance the 
analysis code with die sensitivity analysis capability. 

I. Introduction 

The p- version Galerkin’s method with discontinuity in time developed here was 
based upon the work done by Bey and her colleagues 1 ' 3 . The method uses hierarchical 
high-order polynomials to interpolate the temperature through the thickness of the 
panels. Therefore, the method can catch high thermal gradients without introducing 
geometrical discretization through thickness. This special feature of the method offers 
an order of reduction in finite element modeling. That makes the method very 
attractive to many engineering applications. Furthermore, the method treats time as 
the fourth dimension in its discretization and formulation. The weak continuity 
condition is imposed at the joints between each time intervals (or elements ). The end 
result is a high order, single step time marching scheme. Though the stability of the 
method is ensured, the enror analysis and the adaptivity of the method have been the 
focus of recent research. 

The emphases of the work reported here are placed upon the numerical modeling 
of temperature-dependent and distributive nature of material properties. Such 
approximation not only affects the accuracy of the analysis but also the formation of 
Jacobian matrices which are required in the Newton-Raphson’s method and 
sensitivity analysis. Numerical examples are provided to show the relations between 
the modeling parameters and the accuracy of the method. An attempt is also made to 
analyze a heat conduction problem with a moving heat source. 
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Sensitivity analysis is defined in this report as a process that derives sensitivity 
equations to compute the derivatives of responses or states with respect to specified 
variables. Since the derivative information can greatly enhance the robustness and 
accuracy of curve fitting, sensitivity analysis becomes a necessary element in many 
engineering applications; such as design trade-off, weather prediction, analysis error 
correction, model adjustment, reliability analysis and design optimization. It is shown 
that the sensitivity equation can be derived and solved using the direct differentiation 
method with a minimal effort in the framework of /^-version Galerkin’s method with 
discontinuity in time. The design variables used in this study are related to material 
distribution. The applications of such derivatives are also briefly discussed in the 
report. 

The report is organized as follows. Section II details the ^-version Galerkin’s 
formulation. Section III discusses the construction and the solution strategies of the 
required finite element equations. Section IV derives sensitivity equations based upon 
the direct differentiation strategy. Section V reports the results of numerical studies. 
Concluding remarks are given at the end of the report. Some of the results have been 
published as a proceeding paper in the 14 th Annual Thermal and Fluids Analysis 
Workshop, August 18 th to 22” , 2003 in Hampton, VA. The paper is listed as the first 
attachment of the report. The computer code developed in this study, along with the 
input of a 3D example, is attached in a CD. 


Il.p-Version, Time-Discontinuous Galerkin’s Method 

The governing differential equation of a general heat conduction problem is given 
as 


du 3 3 a 

pc 2 . E 

dt 

u=f(x,t) 


L, 


du 

l] ~d^ 


= Q(x,t), in Qx(0J] 
on x (0,rj 


33 du r i 

1=1 j =1 dxj 


on 5Q x(0,J] 


3 3 du r l 

*, ]=-h(u-T.), 

1=1 j= 1 dx ; 


on 5 Q a x(0,T] 


and 


« = g(x) , 


in Q at t = 0 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


where the temperature, w(x,t), is the only unknown. Equation 1 represents an initial- 
boundary value problem with Eqs.(2-4) being the temperature and the heat flux 
boundary conditions, and the thermal convective condition; respectively, and Eq. (5), 
the initial conditions. It is assumed that the heat source, Q, the prescribed 
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temperature,/ the flux, q St the thermal film coefficient, h, and the initial value, g, are 
all with proper regularity. In case of material nonlinearity, the specific heat, pc{ u), the 
film coefficient, h(u), and the thermal conductivity, &jj(u) are assumed to be functions 
of temperature. 

Let / = 0 on Q u . Otherwise, u in Eqs.(l-5) can be replaced by u -/to achieve a 
homogenous boundary condition on Q u . As a result, the weak form of Eqs. (1-5) can 
be derived, based upon the Galerkin’s Method for an arbitrary function, w(jc), as 

r, du 2, 2, du dw. e 
I (pc — w+ E E k s )dv+ ihuwds 

a & H>1 dxj an. 

= ^Qwdv+ J q s wds+ jhT^wds (6) 

n an, an* 

In the approach of time-discontinuous Galerkin’s method, the time coordinate is 
treated as the same as the spatial coordinate. The time space is also discretized into 
elements or intervals. Focusing on one time interval, /„ = [*„./, t n ), the weak form, Eq. 
(6), can be extended to the entire product domain Oy.I n 

’/'Hr ^ ) l l huwdsJdi 

n - 1 n UA, i UX j «-l Xl h 

n 

= J \[fQwdv + J q s wds+ ^hT x wds]dt (7) 

where, w(x,t) is the testing function. Furthermore, to enforce the continuous 
requirement at the interface of the time interval at t„.i , a weighted integral form of a 
constraint is imposed to Eq. (7) as 

|(pc + u + - pc~u~)w + dv = 0 (8) 

n 

which is a weak form of the jump condition, Pn-i c n-i w n-\ = Pn-i c n-i w n-n at tn-i , where 
«»-i =limu(x,t„_ 1 -e) 

s — >0 

“n-l = limu(x.t n _ x +e) 

£ — ►O 

and other variables with the superscripts “+” and and the subscript have 
similar definitions. Note that p~_ x c;_ x is evaluated at in i n _ x , whereas 
Pn-iK-i and w »-i 316 evaluated at the same time instant, t n _ x , but in I n . 

m. Finite Element Equations and Solution Strategy 
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The domain of an triangular element is defined as x /„ , which Q e can be further 
divided into two spaces; n p (x,y) and £i d (z), where (x, y, z) is within a triangle 

bounded by three vertices (x^y^ixj.yj) ,and (x k ,y k ) and ~<z<^ where d is the 

thickness. The non- dim ensional local coordinates selected to represent (x,y), z, and / 
are defined by the area coordinates (l,-,!,,!*), and non-dimensional variables, £, and 

x respectively, where 

x = L i x i + LjXj +L k x k 
y = L iyi + Ljyj + L k y k 
1 = 1,- +Lj+L k , 

and 



T=- 5^-1 

tn ~t H - 1 

The solution, u„ = u(x, y, z,t) in an element, Q p *n d xl„, is interpolated in terms of 
^i(Li,Lj,L k ),y/(4) and 8 k (z) as 

U n=Oli 4>iVjQk a ij k ( 1 0 ) 

/=1 y=l*=l 


or in a vector form 

u n (x,t) =x T (x,y,z,t)a K 

=(fix,y)® y<z)<g> 6{t)Y a n (1 1) 

where the symbol, ® , represents the outer tensor product operator and the vectors, tf>, 
y/ and 0, represent collections of basis functions. Particularly, the through-thickness 
basis functions, y/, are made of the Legendre polynomials of the first kind 1 , the 
temporal basis functions, 6, the integrations of the same Legendre polynomials and 
the in-plane basis functions, <f>, chosen for triangular elements, as described by 
Reference 4. 

In Eq. (10), I p , I d and I t represent the numbers of shape functions for planar 
triangular, through-thickness and through-time interpolation, respectively. In our 
case, 

I p =(3 + 3 *(I pp -1 )+(I pp -1 )*(I PP -2)/2) 

I,=(I P ,+V 

where I pp ,Ipd and I pt are the orders of hierarchical polynomials used for 

interpolations. Particularly, in the planar space, there are 3 for vertex shape functions, 
3 * (I pp -l) for edge functions and (I pp -l)* (I pp -2)/ 2 for the interior bubble 
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functions. The details of the interpolation functions for temperature are listed in 
Appendices I to III. 

It is assumed that in this study, the film coefficient h is set to zero and the 
relations between the other material properties, pc and K tj , and the temperature are 

defined in a tabulated form, presented as a result of experiments. Therefore, the 
material properties cannot be explicitly specified as functions of position and time as 
required by integration. An approximation is thus introduced to overcome such a 
difficulty. The standard Lagrange polynomials are used here for this purpose. 

The values of the material properties at the Lagrange points are taken from a 
given material table based upon the values of the temperature found at those points. 
The values of the material properties at elsewhere in an element are then obtained 
through interpolation. In this way, the material properties can be explicitly 
approximated as functions of position and time throughout the problem domain. The 
detailed Lagrange polynomials are given in Appendices IV and V for interpolating 
material properties. 

As an example, the material property, say pc, can be interpolated in an element 
(Op* I n ) as 

pc(x,f) = (N c ( x, y, z, t)) T Sc 

-(N^x.y^N^z^N^Y S ( 12 ) 

i j k 

where, 8 cijk is a component of the vector, 8 C , which takes the value of the material 

properties found in the material table based upon the temperature evaluated at the 
point ( Xj,y it Zj,t k ). Specifically, if the point is called point m, the temperature at point 

m, T m ( x i ,y i ,z j ,t k ), is given as T m = xl> a n where «„ is the unknown of the matrix 
equation in time interval, /„, and = tfx; ,y t )® y( zj )®0(t k ) evaluated at point m. 

Since the code studied here considers the material properties are input in a 
tabulated form. And the relationship between the material properties and the 
temperature is expressed by a precisely linear function as shown in Fig. 2.1. 



Figure 2.1. Material Property-Temperature Relation 
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Thus, if the calculated temperature, T m , falls between a pair of given data points, 
Ti<T m < Tj , the material value, S m , which represents fiqjk in Eq. (12), is given as 
S/ -S t 

S m =S i+ (^—^)(T m -Ti) 

1 j 

or rewritten as 


SjfTj -T i )-T i (Sj -Si) S j -Sj 
" Tj-Ti Tj - Ti 

SjTj-TjSj (Sj -S^ r 

Tj-^ (Tj-Ti) m 

which can be abbreviated as 


)T m 


Pm +a m^m 


where 


SiTj-TiSj 

Pm 


a m 


T -T. 

1 J 

bjSi) 

{Tj-Ti) 

Therefore, one has a linear relationship between the material property and the 
temperature, 

5 m =Pm+(<*mXm) T ‘*„ (13) 


The differentiation of S m with respect to a„ is then obtained as a constant vector, 

06 T 
-^ = ( * m Xm )T 


(14) 


The above equation is needed in computing the Jacobian matrices for Newton- 
Raphson’s iteration. 

Once the interpolation functions are selected and the nonlinear material properties 
are approximated as explicit functions of position and time as shown in Eqs. (1 1-12), 
one can proceed to integrate the terms in Eqs.(7-8) to construct the equivalent 
matrices. The resultant finite element matrix equation for the time interval, I„ , can 
then be expressed as 

l C n (<0+ K n («„ )+ («, )K =9„+ Wn-l +h (15) 

where the subscript, n, indicates that the associated matrix or vector is evaluated with 
functions defined in time interval, I n . Note that each term in Eq.(15) is corresponding 
to an integral in Eqs.(7-8), which can be evaluated based upon the interpolation 
functions of Eqs. (1 1-12). As an example, the capacitance matrix is given as 
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c.=Ic 0 , 

p 

=2 

The details of an elemental capacitance matrix, C n ^ and other related elemental 
matrices are given as Eq.(VI.l) in Appendix VI. 

The matrix equation derived above is based upon the discontinuous Galerkin’s 
method and will be solved in a time-marching fashion. In other words, the matrix 
equation of Eq.(15) will be solved for a time interval at a time, with a n as unknown 
and a n .j known quantities. Because of its nonlinear nature, Eq.(15) can be solved by 
the Newton-Raphson’s method, which leads to a recursive formula 


[c, )+ k, )+ m; (a ';‘ ) + j r ) + j t (a 1 ; 1 > + j, (a';‘ nt*‘, 

= -R~ I 


(16) 


and the unknown is updated by 



+ A a‘ 


(17) 


In the above equation, A a n is the improvement of the solution and R' n 1 is the 
residual of the nonlinear equation at the i-1 iteration, which is defined as 


* = [c. («:" )+ K. (U 1 ; 1 ) + M; (a 1 ; 1 (« )]a 


Moreover, the derivative matrices, J c , and J m are obtained by differentiating 
the coefficient matrices C„ , K n and M+ with respect to the unknown vector a„ , 
respectively. This is usually accomplished at the element level. The derivative 
matrices, J c , Jk and J m are detailed in Appendix VI as Eq.(V 1.2). 


The main steps in the computer code developed for solving Eq.(16) can be 
devised into two major steps, as shown in the flow charts, Figs.3.1 and 3.2. The 
integration of all necessary sub-matrices is conducted in the preparation stage, before 
starting the Newton-Raphson iteration. For example, the three integrals defined in 
Eq.(VI.l) are the three dimensional sub-matrices 


^ 00 = l(9*9 r )^xydA 
a p 

dl 1 

B= J (yip 7 )N x dz= )N z d£ 

d\ -i 


(d 2 ~d\) 


n 1 QQ r 

C= \($*0 T )N t dt= ](»•— )N,dr 
n -i dT 


2 
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which can be integrated independently from the solution procedure. On the other 
hand, the elemental matrices as indicated by Eq.(VI.2) are the results of tensor 
product which can be constructed only after the order of approximation and the new 
solution, a„, are known. Thus, the formation of the elemental matrices as well as the 
matrix equation of Eq.(16) can only be done dining the Newton-Raphson’s iteration. 

IV. Sensitivity Analysis 


Since the matrix equation, Eq.(15), is in the form of a static problem, A(x)x=c. 
Differentiation of Eq.(15) with respect to a design variable, b, gives a sensitivity 
equation for a nonlinear transient heat transfer problem as, 

(C+k+m' +j, +j t +j 

db 


f dC_ 

K db 


dK 

+ 

db 


dM + 

db 


■k + 


j 


dM - 
db 


-M 


H—l 


+ M~ 


da n~, , dq 

db db 


(18) 


The derivative matrices appearing on the right-hand side of Eq.(18) are not available 
yet in analysis phase, but needed for sensitivity analysis. Generation of such matrices 
can be tedious and prone to mistakes. A typical derivative matrix, dC I db, is given in 
Appendix VI. Fortunately, because of the nature of approximation, those derivative 
matrices can be obtained in the same way as the original matrices themselves. This 
becomes evident by comparing Eqs.(VLl) and (VI.3) in Appendix VI. Furthermore, it 

da 

is noted that Eq.(18) is a linear equation of — -, which enjoys the same left-hand 

db 

side coefficient matrices as Eq.(16). Thus, the factored matrices saved from analysis 
at the converged stage, can be reused here to solve the sensitivity equation. That 
makes the sensitivity analysis computationally efficient. 

V. Numerical Results 


The numerical experiences on development and application of the /^-version 
Galerkin’s method are reported here in three separated sections. One main focus here 
is studying the accuracy of the proposed method for analyzing an example, nonlinear, 
transient heat conduction problem. The sources of errors may come from the 
approximation in the temperature interpolation, in the nonlinear material relation and 
in the thermal loads, boundary and initial conditions. Section IV. 1 investigates the 
errors due to approximation in temperature and in nonlinearity. Section IV.2 studies 
the heat conduction under a moving point heat source. The heat source is represented 
by a singular delta function. It is expected that the error is resulted from lacking of 
regularity in loading. However, the results summarized in Section IV.2 represent only 
an initial attempt on the subject. Further study is needed in this regard. Section IV.3 
shows the results of sensitivity analysis which demonstrate possible use of sensitivity 
values. 
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V. 1. Study of Accuracy 

Two examples, one is in 2D and the other in 3D, are used to study four aspects 
of accuracy by varying the orders of approximation in time, space and material 
nonlinearity, by increasing the time and by reducing the size of elements. The 2D 
problem is a lxl square with the exact solution being given as 

u e (x,y,t) = xy(x-\){y-\)t 5 (19) 

, whereas the 3D problem is a lxlxl block with the exact solution being given as 

u e (x,y,t) = xy(x-\Xy-\)(2y/ 2 (z)+'Sif/ J (z)+4y/ li (z)+5yf s (z))t 5 (20) 

where y/^z) is an integrated Legendre polynomial of order i, as described in 
Appendix III. The solutions in both cases satisfy the homogeneous boundary 
conditions and zero initial condition. And the material properties, pc and isotropic & , 
used in these cases are assigned by a linear temperature relation, 

pc,k = a+ flu (21) 

where a and fl are assigned constants. 

It should be noted that to achieve the exact solutions as described by Eqs.(19- 
20), the orders of the in-plane temperature, the through-thickness and the time 
approximation, denoted by P xy , P z , and Ptime, should be at least 4, 5 and 5, 
respectively. Since material properties are linearly related to the temperature as given 
by Eq.(21), the orders of the in-plane, the through-thickness and the time 
approximation, denoted by M xy , M z and Mtime, should be at least 4, 5 and 5, 
respectively, as well. The results reported here focus on the accuracy of the solutions 
which are subjected to varied orders of such approximation. 

Reduction of Errors due to Increase of Order of Approximation in Time 

For linear transient heat conduction problems with smooth solutions, the L 2 - 
norm of the error is bounded by 


INI I 2 ("Ox/ J 


<Chd (Ptime * x) 


( 22 ) 


where the pointwise error, e = u - u e , and u and u e denote the numerical and the exact 
solutions, respectively. The L 2 - norm error is defined for a 2D case as 
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fit 2 = 


fjff- 

0 i U. 


" 2 r &) 2 


, 1/2 


+e 


Licrfy 


<* 


(23) 


The corresponding 3D case can be defined in a similar manner. 

The log-log plot of Fig.5. 1 . 1 for the 2D linear transient problem shows that 
the L 2 - norm of error is linearly proportional to the order of approximation in time, 
(Ptime+1), as indicated by Eq.(22). In addition, the results of the 2D nonlinear transient 
problem shown in Fig.5. 1.2 show an identical, linear relation. Similar trends are also 
observed in the 3D linear and nonlinear cases as shown in Figs. 5.1.3 and 5.1.4. Note 
that in the nonlinear cases studied here, p =1 0 in Fig.5. 1 .2 and 0=1 in Fig. 5.1.4. 
Furthermore, the order of approximation in time for material distribution is 
maintained at 5 for both nonlinear cases. 


To investigate the change of errors along with time, the pointwise error, 
H ! ( Q ), of time, defined below for a 2D case, is used hereafter. 


\ e ( t A l p(cii | jj 


n 


de ) 2 

uJ +< 


dxdy 


1/2 


(24) 


The results of the linear case are summarized in Figs. 5.1.5. Figure 5.1.5(a) shows the 
increase of errors for a longer simulation time period. The effects of orders of 
approximation in time on the error distribution within a single time interval are 
detailed in Figs.5. 1.5(b) and (c). Similar trend can be observed in the 3D case as 
shown in Figs. 5.1.6. It is noted from Figs.5. 1.5(a) and 5.1.6(a) that the nonlinear case 
has a bigger jump of temperature discontinuity crossing the interfaces of time 
intervals than the linear case does. 


Reduction of Errors due to Reduction of Element Sizes 

For linear transient heat conduction problems with smooth solutions, the H 1 - 
norm of the error is bounded by 

(25) 


The results of two analyses with different meshes are summarized in Fig.5. 1.7 for 
both linear and nonlinear cases, where the square-root of element area, -jArea , is used 
to represent the element size, h. The analysis domain is meshed into 4-element and 
16-element cases, as shown in Fig.5. 1.8. Another parameter studied here is the orders 
of polynomials for in-plane approximation of temperature as denoted by Pxy for the 
nonlinear case and LP xy for the linear case in Fig.5. 1 .7. The results show that Eq. 

(25) is closely observed for the cases when the orders of polynomials are 2 and 3 for 
both linear and nonlinear cases. Note that the correct solution is achieved, if the order 
of polynomials reaches 4. 




11 


Increase of Errors due to Increase of Nonlinearity 

Two 3D, nonlinear transient problems are solved with f}=\ and p=A which 
represent different degrees of nonlinearity in material-temperature relation defined in 
Eq.(21). Figure 5.1.9 shows that the degree of nonlinearity can aggravate the error 
caused by inadequate approximation. In this study, the only parameter varied is the 
order of polynomials used in time domain approximation. 

Increase of Error due to Increase of Simulation Time Period 

In this study, the time step is set to be 1 second and the total operational time, t, is 
set to be 4 seconds. The orders of in-plane and the temporal polynomials are selected 
to be 4 and 6, respectively, for temperature interpolation as described by Eq.(ll), 
whereas the orders of in-plane and the temporal polynomials are 3 and 5, 
respectively, for material property interpolation, Eq.(12). The total Lagrange points 
for the in-plane material interpolation is 10 as marked in Fig.5.1.10. Moreover, the 
error is measured by the Z, 2 -norm as defined by Eq.(23). 

At the end of the first time interval, /, = 1 second, the error is in the order of 10" 4 . 
The error is growing with the time. At t 2 = 2 second, the error grows to the order of 
10" 1 . At t A = 4, the error becomes 13.31. The major source of the error is expected to 
be the deficiency in the in-plane interpolation of materials. Since the material 
property is linear in temperature, it’s exact order of interpolation has to be 4, equal to 
the order of temperature interpolation. However, in this study, the order used is 3. 
Such kind of error can be aggravated from one time interval to the next. The results of 
the exact temperature distribution are shown in Fig. 5.1.11. 

V. 2. Heat Transfer Analysis with a Moving Heat Source 

An attempt is made here to simulate the welding process in which the testing 
article is subjected to a moving point heat source. Let the point heat source move 
along with y-axis with constant velocity, v, and an intensity, q. Thus, the point heat 
source is represented as 


Q = q*8(x-x,y-y) (26) 

where in our study, x is fixed and y =tv. The equivalent heat source vector can be 
obtained by the following equation, 

f = J \q*8(x-x,y-y)xdAdt 

In&e 

-JJ* x8(x-x,y~ y)[<p(x,y)<gi0(t)]fIAdt 

= J q*]p(x,y)®0(t)\it 

In 


( 27 ) 
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The area coordinates corresponding to x and y can be found in a triangular 
element defined by three vertices, (jc„ y,), i=\ to 3, as 

x = L l x 1 +L 2 x 2 + Z.3JC3 
y = tv = L l y i +L 2 y 2 +L 3 y 3 
1 = Zq + lt 2 + JL 3 

with which f(x,y) can be expressed in terms of t explicitly. Equation 27 can then be 
integrated to obtain the equivalent thermal load vector in a given time interval. 

Two meshes are studied here; 4-element and 48-element models as shown in 
Fig.5.2.1. The temperature distributions resulted from linear simulation are shown in 
Figs.5.2.2 and 5.2.3 at time instances equal to 0.1 second and 1 second, respectively, 
based upon a 4-element mesh. The advance of the high temperature zone is evident 
from these two figures. However, the movement of the highest temperature spot is 
about 4 time faster than that of the point heat source. Furthermore, the highest 
temperature obtained from a 48-element mesh, as shown in Fig.5.2.4, is about 40 
percent higher than that obtained from a 4-element mesh. These observations make 
the simulated results unacceptable. It is expected that the singularity of the point load 
given by Eq.(26) causes numerical difficulty. Further study is needed to gain better 
understanding of such type of application of the /^-version Galerkin’s method. 

V. 3. Sensitivity Analysis 

The material properties studied here are considered temperature-dependent. 
Therefore, any parameters that define the material-temperature relation can be 
considered as a design variable. In Case 1, the slope of such relation is considered as 
a design variable. Furthermore, since the temperature is modeled as a distributed 
function within an element, the material properties here can also be viewed as 
distributive. Therefore, the material distribution within an element can be viewed as 
a design variable as well, regardless whether the problem is nonlinear or not. In Case 
2, the material properties - die thermal conductivities at the Lagrange points of the 
entire problem are considered as design variables. In order words, it is the material 
distribution considered as a design variable in Case 2. 

Case 1 


The slope of the thermal conductivity-temperature relation is considered as a 
design variable. Since only the matrix K depends upon this special design variable, 
the right-hand side of Eq.(18) can be greatly reduced to a single term, (dK / db) a. The 
results of the thermal sensitivity coefficients at times from 1 second to 4 seconds are 
shown in Fig.5.3.1. Note that the thermal sensitivity coefficients is interpolated in the 
same way as the temperature; i.e., 


du 

~db 


= X T (x,t) 


da 

lb 
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where da / db is obtained from Eq.(18). Comparing with the finite differencing, the 
errors of the thermal sensitivity coefficients are less than 1(T* for all the time 
intervals. 

Case 2 


Since in this study, the material property is assumed to be a distributed function, 
one may then assume that the square slab is made of various materials. In this 
particular case, the value of k at each Lagrange point is determined by it’s own 
material table. If the slope of each material table is considered as a design variable, 
there are 10 independent design variables in total, as marked in Fig.5.1.10. 

Figures 5.3.2 and 5.3.3 show the distributions of the thermal sensitivity 
coefficients of du / dbi and du / dbg , where bj and bg are the slopes of the thermal 
conductivity-temperature relations at Lagrange points 1 and 9, respectively. The 
figures reveal that the design variable, bj t affects the change of temperature along the 
diagonal line, whereas the design variable, bg, does the same to the area off the 
diagonal line more. Finally, all 10 thermal sensitivity coefficients of the temperature 
at the center point are collected and plotted out in Fig.5.3.4. The picture indicates the 
degree of influence of individual design variable on the temperature at the center 
location at different time instances; i.e., du(x c ,y c )/ db t ,i = 1 to 1 0. It is as expected that 
the figure shows that the value of the thermal conductivity at the center affects the 
temperature at the same point the most. 

VI. Discussions and Remarks 

The report documents the effort in the development of a structural-compatible 
code for nonlinear heat transfer analysis and sensitivity analysis based upon the 
framework of the p- version time-discontinuous Galerkin’s method. Numerical studies 
are made to assess the performance of the method. 

It is found that the error characteristics of the method in solving the nonlinear 
problem are very similar to that of the linear one, if the degree of nonlinearity is 
small. The high degree of material nonlinearity can aggravate the numerical errors 
through the approximation in material distribution. Furthermore, a better approach is 
needed to convert general, pointwise descriptions of the loading, boundary and initial 
conditions into equivalent and consistent vector forms. These concerns point to the 
need of developing an adaptive method for error control. Note that the leading 
elemental matrix in the transient problem is non-singular. This observation makes the 
element-by-element error correction scheme 5 developed for linear, static problems 
very attractive to the nonlinear, transient problems studied here. 

An attempt is also made to find the thermal sensitivity based upon the /^-version 
time-discontinuous Galerkin’s method. The parameters that are related to the 
material-temperature dependent relation as well as the material distribution are 


14 


considered as design variables in this study. Though construction of the matrix 
equation for thermal analysis is complicated, construction of that for thermal 
sensitivity analysis is rather simple. Solving the resultant sensitivity equation is also 
demonstrated to be computationally efficient. 

Although the p-version Galerkin’ s method can effectively achieve an order of 
reduction in meshing effort, the size of the matrix equations remains large for solving 
a 3D problem. An effective iterative solution strategy for a non-symmetric matrix 
equation is needed in order to support more challenging applications of the proposed 
method. 
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Fori to# of Time 
Elements 


Figure 3.1. The Prepration Stage and Marching through Time Intervals 
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End Newton-Rapson 
Iterations 


For 1 to # of 

Newton-Rapson 

Iterations 


Figure 3.2. The Newton-Rapson Iteration for a Single Time Interval 











Number of Time Steps 


Figure 5.1.1. Error, ||ej| L2 , versus Number of Time Steps; 
2D, Linear Case, P = 4 at t = 2 seconds 
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Figure 5.1.2. Error, ||e|L .versus Numberof TimeSteps; 

2D, Nonlinear Case,P xy = M xy = 4 and P = 10att = 2 seconds 
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-S' P tiTB =5 


Nurter of 11 me Steps 

Figure5.1.4. Error||e|| L2 versus Time; 3D, Nonlinear Case, 
P xy = M xy = 4, P z = 6, M z = 5 at t = 2 seconds 
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fel sac t=2 sac t=3 aac t=4 sec 

Figure 5.1.1 1. Exact Temperature Distribution 





Figure 5.2.2. Temperature Distribution at t=0.1 Second for a 4-Element Model 
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Figure 5.2.3. Temperature Distribution at t=l Second for a 4-Element Model 
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0 0.00190909 0.0320626 0 0890424 0.179123 0.32538 



t=i sec t=2 sec t=3 sec t=4 sec 


Figure 5.3.1. Coefficients of Thermal Sensitivity with Respect to Homogeneous 

Thermal Conductivity; . 

db 
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0 2 .82828E-050 .000532121 0.00169899 0.00438586 0.00772727 



Ume= 1 tim»=2 timer 3 *n#=4 


Figure 5.3.2. Coefficients of Thermal Design Sensitivity with Respect to Thermal 

'V / v 

Conductivity at Material Lagrange Point 1; — — 1 — - . 

db , 


38 



0 0.000247475 0.00476061 00169394 0 0515859 0 100152 



1=1 9BC t=2 HC t=3 MC t=4 flee 


Figure 5.3.3. Coefficients of Thermal Design Sensitivity with Respect to Thermal 
Conductivity at Material Lagrange Point 9; . 
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Appendix I 

Polynomials for Temperature Interpolation (up to 6 order) in XY Plane 

T=(p®0®y) T a 


Comer Functions 


4=^1 
^2 =I 2 


Edge Functions: (Pi p -1) per edge 


h'S-Lj-L, 

<j> 2e = ~v/l0 • {Lj ■ L 2 k - Lj ■ ) 


fae ~ '($ A/ '4 10-Z.y -l} K +S-I?j -L k — Lj -L k ) 

(f» Ae = — ■ ( — 7 Lj -L k + 21 • l}j ■ L k -2\- L 3 j -I? k +7 • L*j ■ L K + 3 • Lj ■ l} K - 3 -L 2 j -L K ^j 

<tse = -21 -Lj L s k +84-4 4 -126-4 4 +84-4 4 -21*4 Ac 

+ 14 I j - 4-28-4'4 +14-4-I Jf -Z y -Z, jr ) 


e 


k 

1 

1 

2 

2 

2 

3 

3 

3 

1 


Permutation table for edge number, e, and comer points, j and k. 




Bubble Functions: (Pi p -l)*(Pip-2)/2 


0\ b — A ■ 
02 B ~ L\ 
03B ~ A 

04 B 

0SB ~L\ 
06B ~ L\ 
01 B = A 

08 B ~ A 

09 B = A 

= A 


L^Lj, 

• £ 2 • Lj • (Lj - ij ) 

■ L 2 ■ L 3 -(p.- L$ — l) 
L 2 L,(L 2 -L l )-(2-L 3 -l) 

■L 2 -L i y^-(L 2 -L 1 f-l) 

L 2 -L 3 j^(2-L 3 -lf-l) 

^L 3 ~^{L 2 -L i y-l){2L 3 -l) 

■L 2 .L 3 .(L 2 -L,)-^-{2L 3 -lf-l) 

-L 2 .L 3 ±-I?-(L 2 -L i y-3.(L 2 -L l j) 

• L 2 -L 3 -y(s(2-l 3 -lf~3(2-L 3 -l .)) 
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Appendix II 

Lengendre Polynomials for Temperature Interpolation (up to 6 order) 

in Time 


Lengendre Polynomials (- 1 < x < l) 

e x =t 

e,.i(sr’-3r) 

0 A =^35r 4 -30r 2 + 2) 

O 

0 5 =^63r 5 -70r 3 +15r> 

O 

0 6 =—(23\ t 6 -315r 4 +105r 2 -5; 
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Appendix III 

Polynomials for Temperature Interpolation (up to 6 order) 
through Thickness in Z 

Integration of Legendre Polynomials (-7 < z < l) 


Vo 




Vi 

-i a ^> 




1 ,3 2 

3 . 


V 2 

‘Vs'V* ' 



v 3 

-4-< 5 -s- 

Vio 2 

i z) 



1 .35 4 

42 2 

7 . 

V 4 

= —==( — z’ 

Vl4 8 

z 

8 

+ 8 ^ 


1 .63 5 

90 3 

27 . 

Vs 

= ~r=( Z 3 

■M 8 

z 

8 

H z) 

8 


1 .231 

6 385 

4 165 

Vs 

= —f=( 2 

V 22 16 

16 

z + 

16 
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Appendix IV 

Lagrange’s Polynomials for 
Interpolation of Material Properties in XY Plane 


Order 1, 



Comers 

tf,i = A 
N p 2 =L 2 
N p 3 =L 3 


Order 2, 



Comers 

N P i = ( 2L \ ~V'h 

N pl =(2L 2 -l)L 2 

N p3 =(2L 3 -\)-L 3 

Mid-points 

N pi =4-L l L 2 
N p s —4- L 2 L 3 
N p 6 = 4 • 



Order 3, 



Comers 

N pl =\(3L l -l)-(3L l -2).L l 

N p2 =±(3L 2 -l)(3-L 2 -2)L 2 

^3=|(3i 3 -l)(3I 3 -2)^ 3 

Mid-points 

* P 4=|-W3-A-1> 

^5=|-i 1 L 2 r3-i 2 -i; 

N p6 = -L 2 L 2 (3-L 2 -l) 

= — •L 2 L 3 (3-L 3 -1 ) 
N pt ^~L 2 L l (3L 2 -\) 
N p9 =~L 2 L x (3L l -\) 

Center 


N p iq — 27 ■ L^L 2 L 2 
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Order 4, 


3 



Comers 

N p i=^(4L i -\)(2L 3 -\)(4L l -3)L y 
X P 2 =\(*' L 2 -V-(2 L 2 -\)(4 -L 2 -3 ) L 2 
N p i = ^(4 • i 3 - 1)- (2 ■ L 3 - 2> (4 ■ is - 3) • L 3 


Mid-points 

^=“W4-i,-Ur2-i,-i; 

N p5 =4L l L 2 (4L l -l).(4-L 2 -l) 
N p6 = --L l L 2 (4-L 2 — 1 )-(2-L 2 — 1 ) 

N p1 = jL 2 L 3 (4L 2 -1)-(2L 2 -\) 
N pt = 4 ■ L 2 L 3 (4 ■ L 2 -\)-(4-L 3 — 1 ) 

N p9 -~L 2 L 3 {4-L 3 -l) (2 L 3 — 1 ) 

Hpto = ~^'^ 3 L l (4- L 3 -1 )-(2- L 3 —1) 
Npu =4 L 3 L 1 (4 L 3 - \)-(4-L l -1) 
N p n =— L 3 L x {^ L i -1 ) (2 L l - 1 ) 

Center 

N pl3 =32-L l L 2 L 3 (4-L,-l) 

N pU -32- L l L 2 L 3 (4- L 2 -l) 

^ pis —22-L^L 2 L 3 (4-L 3 —\) 
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Appendix V 


Lagrange’s Polynomials for Interpolation of Material Properties in Time 

and Through-Thickness 


1 st order 


N(£)=tf 0 (<f)-K 0 


2 nd order 


W 0 (£)=~(W) 


N(?)=tf 0 fc).« 0 +*,(£).«, + N 2 {g)u 2 



1+S) 

n 2 {§)= W 2 

3 rd order 


^)=N^)-u 0 + N^)- Ul +N 2 {s)-u 2 +N 3 {z)-u 3 


-1 -1/3 1/3 1 

N 0 (^)=~{-l + 4 + 9^ 2 -9-^) 

10 

^ife)»77-( 1+^-9 ^ 2 -9^ 3 ) 

10 

N 2 ^) = ~{9-21^-9-e ^21-e) 
10 

N 3 {/;) = ~^ + 21-l;-9-Z 2 -21-?) 
10 
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4 th order 


Nfe)= +AT 1 (f)-* 1 +N 2 (Z)-u 2 + nM» 3 +N 4 (4)-u 4 


-1 - 1/2 0 1/2 1 


^ o (^)=~(#-^ 2 - 4 -^+ 4 -^) 
o 

V&h—fe+Z 3 - 4 ^ 3 " 4 ^ 4 ) 

y -( 8 -^- 16^ 2 - 4-^ 3 + 16 -<? 4 ) 

JV 3 (£) = i -(8 ^ + 16^ 2 - 4-<? 3 - 16 -<* 4 ) 
^ 4 (^) = (l-5-# 2 + 4^ 4 ) 

5 th order 


NfeMofeK +*i(f)-«. +^ 2 fe)-«2 + *,($)•«, +tf 4 fc)-«4 +^ 5 fe)-«5 


-1 -3/5 -1/5 1/5 3/5 1 


^ o (^) = — ( 9 - 9 -^- 250^ 2 + 250 -£ 3 + 625-£ 4 - 625 -£ 5 ) 

768 

= — ( 9 + 9 -£- 250-£ 2 - 250-£ 3 + 625-£ 4 + 625 -| 5 ) 

768 

^ 2 ( 1 )= — ' ( 75 - 125 - 4 : - 1950 - 4 :2 +- 3250 -^ 3 +1875 <^ 4 - 3125 -^ 5 ) 
768 

W 3 fe )=— • -( 75 + 125 -£-1950 £ 2 - 3250-£ 3 +1875 £ 4 + 3125 -£ 5 ) 
768 

AT 4 (^)=— -( 450 - 2250 -^- 1700 -^ 2 + 8500-£ 3 + 1250-£ 4 -6250 £ 5 ) 
768 

N 5 ( t )=- L .( 450 + 2250 -^ - 1700-£ 2 -8500 £ 3 + 1250-£ 4 + 6250 -# 5 ) 
768 
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Appendix VI 

Sample Matrices for Analysis and Sensitivity Analysis 


* 

*. 2 


c„. = I I J pcX^-dzdAdl 

2 

d 

2 ( d(t\ T 

= J J f — dzdAdt 

t.-, o t zL V dt ) 

2 

d 

2 *. AO T 

= 111 ^ WN^dA® l^N^dz® }e— N^dt)dy k } 


i J * a , 


-i 

2 




dt 




Vitiprl 
2 

d 

». 2 


2 




adzdAdt 


d0 T 


= ZZZ^J’ > ’ >Tjff c P i dA ® Jf^czkdz® j®— W^dtla ](a mXm )} 

* J k Hp 


-d 


tn-j 


dC, 


db 


d_ 

2 

u ni 

", z' 


db St 


dzdAdt 


= '] J )(N T e ^X<l>®V®9{<l>®V ®— 1 <fedL4<ft 

2 


J# r AVt4® }rr r ^«.*® J o^- N«*)^rr) 


i J * fl. 




* 




(VLl) 


(VI-2) 


(VL3) 



50 


Attachment I 

Paper Presented at the 14 th Annual Thermal and Fluids Analysis 

Workshop 

August 18 th -22 nd , 2003 
Hampton, Virginia 



SENSITIVITY EQUATION DERIVATION FOR TRANSIENT HEAT 
TRANSFER PROBLEMS 


By 

Gene Hou and Ta-Cheng Chien 

Department of Mechanical Engineering 
Old Dominion University 
Norfolk, VA 23529-0247 

Jeenson Sheen 
Department of Technology 
Norfolk State University 
Norfolk, VA 23504 


ABSTRCT 

The focus of the paper is on the derivation of sensitivity equations for transient 
heat transfer problems modeled by different discretization processes. Two examples 
will be used in this study to facilitate the discussion. The first example is a coupled, 
transient heat transfer problem that simulates the press molding process in fabrication 
of composite laminates. These state equations are discretized into standard ^-version 
finite elements and solved by a multiple step, predictor-corrector scheme. The 
sensitivity analysis results based upon the direct and adjoint variable approaches will 
be presented. The second example is a nonlinear transient heat transfer problem 
solved by a p - version time-discontinuous Galerkin’s Method. The resulting matrix 
equation of the state equation is simply in the form of Ax = b , representing a single 
step, time marching scheme. A direct differentiation approach will be used to 
compute the thermal sensitivities of a sample 2D problem. 


INTRODUCTION 

Sensitivity analysis is defined in this paper as a process that derives sensitivity 
equations to compute the derivatives of responses or states with respect to specified 
variables. Since the derivative information can greatly enhance the robustness and 
accuracy of curve fitting, sensitivity analysis becomes a necessary element in many 
engineering applications. Examples include design trade-off, weather prediction, 
analysis error correction, model adjustment, reliability analysis and design 
optimization. 

In sensitivity analysis, the response or dependent variable can be a real number, a 
function or a fimctional and the independent variable can be a real number or a 
function. The challenge of sensitivity analysis arises when the responses to be 
differentiated involve the solutions of some governing equations. In these cases, the 



responses or states are implicitly related to the independent variables through the 
governing equations. These governing state equations can be expressed as 
differential, integral or algebraic forms. The latter is usually as a result of numerical 
discretization of the former. 

Many approaches are available for sensitivity analysis, including automatic 
differentiation 1 , complex variable method 2 , finite differencing, or more traditional 
analytical approaches 3-6 . The analytical approaches used for sensitivity analysis can 
be further classified into various categories. It can be classified as the discrete 
approach or the distributed (continuous) approach, based upon which types of 
governing equations, responses and independent variables are considered. The 
discrete approach works with discretized algebraic equations and real numbers, while 
the distributed approach works with functionals and functions. Sensitivity analysis 
can also be classified as the direct differentiation approach or the adjoint variable 
approach, based upon whether the derivatives of the state variables are computed 
explicitly in the process. The bulk of the effort of the direct differentiation method is 
to establish a sensitivity equation solved for the derivatives of state variables, whereas 
the effort of the adjoint approach is to form and solve an adjoint equation and 
eventually, compute the derivatives of the responses without explicitly computing the 
derivatives of the state variables. 

Subjects related to sensitivity analysis of dynamic transient problems can be 
found in the literature 3 "*' 6 . Thermal transient analysis can be performed in the similar 
fashion 7 " 9 . Recent development and applications of thermal transient analysis can be 
found in aerospace applications 10-11 , laser surface treatment 12 , material processes 1314 , 
and their references. The major focus of the report is on the derivation of sensitivity 
equations for transient heat transfer problems represented in various forms of 
governing equations. Two examples will be used in this study to facilitate the 
discussion. 

The first example is a coupled, transient heat transfer problem that simulates the 
press molding process in fabrication of composite laminates. The state equations are 
made of a heat conduction equation that calculates the through-thickness temperature 
distribution and an empirical equation that monitors the chemical-kinetic reaction of 
resins. These state equations are discretized into standard ft-version finite elements 

and expressed in the form of Ax+ Bx = c . The resultant equations are then solved by 
a multiple step, predictor-corrector scheme. Both of the direct and adjoint variable 
approaches will be used to derive the sensitivity equations in continuous forms. The 
numerical implementation aspects of those sensitivity equations and their accuracies 
will be studied in this paper. 

The second example is a nonlinear transient heat transfer problem solved by a p- 
version, time-discontinuous Galerkin’s Method 15 ' 17 . The resulting matrix equation of 
the state equation is simply in the form of Ax = b, representing a single step, time 
marching scheme. A direct differentiation approach is used to compute the thermal 



sensitivities of a simple 2 D heat conduction problem. Here, possible usages of the 
sensitivity results are presented. 


EXAMPLE 1: A- VERSION FINITE ELEMENT APPROXIMATION 

The interest of this study is in compression molding of a filled polyester resin 
reinforced by chopped glass fibers. The unmolded composite is produced in sheets 
which are from 3 to 6 mm thick, typically. The resin consists of a thick dough and the 
chopped fibers (about 25 mm long) which are randomly oriented in the plane of the 
sheet. In this form, the material is called sheet molding compound, or simply SMC 18 . 
A diffusion reaction system in terms of the temperature distribution and the degree of 
cure can be used to describe the cure process of the composite material under 
consideration. 


The system equations can be expressed as 


du d 2 u da 

pc~ — k— T = pH r —, 
dt dx 2 dt 

in (0,/t)x(0,r] 

(1) 

u{h,t) = u c (t) 

on (0,r] 

(2) 

du(0,t) 

dx 

on (0, t ] 

(3) 

and 



“ = “oM» 

in (0 ,h) at/=0 

(4) 

where p and c are the density and the specific heat of the composite 

material. 


respectively, k is the thermal conductivity in the direction perpendicular to the plane 
of composite material, 2 h is the total thickness, u c (t) is the cure cycle in K ° , T is the 
time in seconds needed for the completion of one cure cycle, H r is the total or 
ultimate heat of reaction, and the last term in Eq. 1, p H r ^ a /^ t , is the rate of heat 
generated by the chemical reaction as characterized by the degree of cure a . 

The degree of cure, a , is defined as the fraction of heat, //(f), released up to time, 
f, for the resin system under cure; a = H(t )/ H r Both H(t) and H r in Eq. 4 can be 
measured experimentally by Differential Scanning Calorimetry (DSC). For an 
uncured material, a approaches zero, and for a completely cured material, a 

approaches one. The reaction rate, . depends strongly on the curing 

temperature. As an example, the cure rate equation of a stepwise isothermal curing 
process which can be used for a polyester SMC is described as follows: 


IT -t<« T > 

dt 

= ( Kj+K 2 a n Xl-a) n (5) 

= ( a 2 e~ d ' ,Kr + a 2 e~ dl /RT a m ){l-af 


where a x ,a 2 ,d l ,d 2 ,m and n are constants, R is the universal gas constant, and K x 
and K 2 are exponential functions of the temperature. 

The optimal cure cycle design 19 aims to select the profile of cure temperature, 
u c (t ) , to achieve the following goals in a compress molding process. 

a) The maximum temperature inside the composite during the cure process can 
not be too high to avoid buring. 

b) The material is cured completely at the end of the cure process. 

c) The material is cured uniformly at any time during the cure process. 

The first two objectives may be mathematically formulated as point-wise functions, 

T(t)<T f in [0,h]x[0,r] 

a(x,T)>a f in [0,/t] 

where x is the total time required to complete one cure cycle. Furthermore, 
the last objective may be measured by the temperature uniformity, which is expected 
to lead to the uniformity of the curing reaction. Mathematically, the temperature 
uniformity is represented by the least-square integral of the deviation between the 
point-wise temperature and the averaged temperature as 


y/ 0 = \[ ju 2 dx-( Judx) 2 / h]dt 
0 0 0 


( 6 ) 


To support the optimal cure cycle design, the thermal design derivatives of the 
functional, y / 0 , the temperature, u(x,t) and the degree of cure, oc(x,t), with respect to 
the cure temperature, u c (t ) , are required. 

Note that the cure temperature appears as a part of the non-homogeneous 
boundary equation in Eq. 2. By using the following replacement of the temperature, 

u(x, t), by u (x, t) as 

u(x,t) = u(x,t)+u c (t) 

which leads to simplification of the heat conduction equation, Eq. 1, 
du , d 2 u du ( - \ 

**= k w- pc ir +pH ' f(a - u+u ' ) 


( 7 ) 


with the homogeneous boundary equations 


u{h,t)~0, 

dx 


in [0,7 ] (8) 

in [0, r ] (9) 


and the initial condition. 


u(x,0) = u 0 (x)-uJ0), in [0,h] (10) 

where /is defined in Eq. 5. Since the initial temperature, u 0 (t ) is the same as the 
initial cure temperature for most applications, Eq. 10 may yield a homogeneous initial 
condition as well. From here on, u(x,t) is abbreviated as u(x,t) for simplification. 


The weak forms of the above equation can now be derived based upon the 
Galerkin’s method for arbitrary functions, w(x) and s(x), as 

hf 3.. ni.. j.. "\ 


n. 


-•-fl 


du d 2 u du 

pc- l c -— + p c £ 

dt Bx 2 dt 


— pHrfydz 


(ID 


and 

^ f)/Y 

n a =0=\(~f)sdx (12) 

o ot 

THERMAL SENSITIVITY ANALYSIS 

The system equation of Eqs. 1-5 simply reveals the fact that their solutions, u(jc,r) 
and a (x,t) are functions of the design variable, u c (t). However, since the design 

variable u c (t) itself is a function, the thermal derivative of u(x,t) with respect to u c (t) 
can be defined as the variation of u(x,t), du , due to the variation in u c (t), 8u c , 

du = ~u( x,t;u c + edu c 

The cure derivative, 6a can be defined similarly. 

The variation of the functional defined in Eq. 6 is then given as 

dyf = \[ j2u - — fudx Jdudxdt 
oo ho 


( 13 ) 



It is assumed that u(x,t;u c ) and a(x,t;u c ) have enough regularity in the time- 
spatial domain and in the design space. Thus, the order of the differentiation and the 
variation is exchangeable; therefore, 

d(&) _ 
dt larj 
d(Sa) _ / 3a N 

dt ^ dt j 
and 


a® 

dx 2 


= J[v* 


With the aid of the above equations, the variations of the state equations defined by 
Eq. 11-12 yield 


T h 


0=11 


00 




A- pH r -^—SaA 
dx 


da 


- pH r f - SuA - pH r 5u c A 

du du , 


dxdt. 


OllU 

0 = ))[sS^^--s^-Sa-s^-du-s^5u\xdt. 
ooV dr da 3u du c J 

Integrating by parts simplify the above equation as 

VJ7 dA ,d 2 A „ ,df\ f dA „ . 


du 


du. 


c / 


-pH r A^-6a 

da 


dxdt+ f k~du dt 
„ J dx „ 


- JfeA d(*^) + dlx-t- Jpc/ldu ( .| dx, 

n dx , n n n n 


(14) 


(15) 


( 16 ) 


and 


( 17 ) 



THE ADJOINT VARIABLE APPROACH 

Adding Eqs. 13, 16, and 17 up, one has the variation of the functional \j/ as 


x Vrl f , 2 \ A 3/1 , d 2 A „ df df 




oo\ 

rh( 




0 0 
h 


dA 


df 3/ 




Su^dxdt 


cj 


T * 22 * r 2/ii 

+ \pcASu c | + Jjfe— — Su\ dt- \kA— — 
o t o OX 0 0 ox 


dt 


Sudxdt 


h * h r 

+ \pcA5u c \ dx + jsSa\ dx 
o ooo 


(18) 


Note that A(x,t) and s(x,t) are arbitrary functions, and the design derivatives Su and 
Sa are the only two unknowns in the above equation. One may now specify the 
variables X and s in such a way that all of the terms associated with Su and Sa are 
dropped. This can be accomplished by introducing the following adjoint equations for 
X and s as: 


n i^A u 

0 = pc— + k—+pH r 


dA 

dt 


3 / Bf 

du du 




(19) 


and 


n ds .df df 

0= ^: +f>H ^T~ +s T- 
dt da da 


( 20 ) 


with the terminal conditions. 


A( 0 ,t) = 0 , 


in [0, h] 


( 21 ) 


( 22 ) 


s(x,t) = 0, 


and the boundary conditions, 


cU 

dx 


(o,/)=o. 


in [0, h\ 

in [0,7 ] (23) 


/l M = o, 


in [0,7] (24) 


Thus, the combination of Eqs, 18-24 provides a simple formula for the design 
derivative of the functional. 


T h i 


dA 


a/ 








| /Ot/l(o)^u c (o)dx 


(25) 


o 

Equation 25 shows that the design derivative of iff , namely, Sl/f , is a functional 
of state variables a and u, and the adjoint variables A and s. The matrix equations 
which solve the nodal values of a and u can be formed by the standard h - version 
finite element method as. 


Ca+ Ka ~ q 


(26) 


Ma = r 


(27) 


Since the adjoint variables of Eqs. 19 - 20 form an “adjoint” diffusion-reaction 
system similar to the original one, the same numerical scheme used to solved the state 
variables a and u can be extended here to compute the adjoint variables A and s. For 
instance, using the same shape functions of u and a to interpolate the adjoint 
variables X and s obtains the following matrix equation for nodal value of A and s in 
the form of 

CA-KA = q * (28) 


and 

Ms = r* 


(29) 


with proper boundary and terminal conditions. 

In general, the adjoint equation cannot be solved simultaneously alongside with 
the original system equation. Because of the terminal conditions, the adjoint 
equations can be solved by either the backward integration along the real time t - axis 
directly or the forward integration along the artificial time r* - axis, provided that the 


independent variable t is changed to t* ax t* = r-t . However, both approaches 
require the solution of the original system equation known prior to solving the adjoint 
equations. 


THE DIRECT DIFFERENTIATION METHOD 

The direct differentiation method is an approach that differentiates the governing 
equations with respect to a design variable directly. The variation Su in <fy/of Eq. 13 
and 5a can be obtained by solving the equations, 8n a =0 and 8n u = 0, in then- 
weak forms such as Eqs. 14-15, or in their differential forms as 


pc 




dt 


dx 


dt 


da 


+ pH T $-8u + pH r $-8u c 
du du „ 


(30) 


and 


dt da du du„ 


(31) 


where the design variable variation, 8u c , is known. 

Again, Eqs. 14-15 or Eqs. 30-31 can be discretized into finite element matrix 
equations, similar to Eqs. 26-27 as Su and da can be interpolated by the same shape 
functions as u and a . These equations can be symbolically written as 


Cab+Ka b =q b 


(32) 


M a u 


= r b 


(33) 


where 8u and 8a are interpolated by the same shape functions as u and a , in which 
a b and a b are the nodal vectors of Su and 8a . 

NUMERICAL RESULTS 

The initial-value problems of the state variables, their thermal derivatives and 
corresponding adjoint variables are all solved by the computer code, DE 20 . The DE 
program is one of predictor-corrector integration algorithms using Adams family of 


formulas. The truncation error is controlled by varying both the step size and the 
order of the polynomial approximation. The DE program is quite easy to use and has 
the capability to manage moderate stiff equations which happen commonly in the 
problems of chemical kinetics. To maintain a unified accuracy in the analysis, the 
computation of two state variables, namely, the temperature and the degree of cure, 
are subjected to the same error tolerance in this study. 

Note that the coefficient matrices of a b and a h in Eqs. 32-33 are identical to 
those of a and a defined in Eqs. 26-27. Therefore, the same numerical scheme and 
the numerical tolerance can be applied to solve both Eq. 26-27, and 32-33 
simultaneously for state variables, a and a , and design derivatives, a b and a b . In 

this way, the state variables and the design derivatives can enjoy the same numerical 
accuracy. 

Regarding the computational efficiency, it is worthwhile mentioning two notes 
here: 

a) Because the coefficient matrices of Eqs. 32-33 are identical to those of Eqs. 26-27, 
the triangular factorizations of matrices C and M are needed to be done only once. 

b) Compared to the original system equations, the right hand sides of the equations 
for computing a b and a b , such as Eqs. 32-33, may have different frequency contents. 
Thus, to maintain the same numerical accuracy, a small time step At may be required 
for the DE program to solve the pairs (a,a ) and (a b ,a b ), simultaneously. 

Once a and a b are available, the design derivative given in Eq. 13 can be easily 
obtained by numerical integration. Another suggestion is to rewrite the integral form 
of Eq. 13 as a differential equation of Slff given as 


dSy/ 

dt 


= J| 

« o 


Sudz 


(34) 


The above ordinary differential equation of Sy/ b can then be solved simultaneously 
with equations of {a,a ) and (a b ,a b ). In this way, one extra equation of design 
derivative Slff for each design variable is added in the design sensitivity analysis. 
However, the accuracy of Sy/ is secured. Eq. 34 is used to generate the current 
numerical results. 

An example which deals with the cure process of compression molding is 
presented in this section to discuss the numerical accuracy of the adjoint variable 
method and the direct differentiation method for calculating the thermal design 
derivatives. The accuracy of the thermal design sensitivity can be checked by using 
the fundamental definition of design derivatives which can be approximated by the 
finite difference. 


The finite perturbation of the design variable, Ab , is defined as the difference 
between a perturbed design b’ and the nominal design b, i.e., Ab = b* - b. As a 
result, it follows that 

Ay/ = w(b')-w{b) (35) 

= yr Ab 


The above equation provides a simple means to check the accuracy of the design 
sensitivity analysis. Nevertheless, the difficulty of this method is the selection of Ab . 
If Ab is too large, the approximation in Eq. 35 is not valid. On the other hand, if Ab 
is too small, the round-off error in the computation of [y/(b ’ ) - y/{b)\ becomes too 
large to ensure the validity of Eq. 35. 

The first example present here deals with the cure process in which the cure 
temperature of the process is assumed to be a constant temperature. The nominal cure 
temperature is taken as 423° k . According to the approximation defined in Eq. 35, the 
results oi Ay/ and Ab shown in Fig. 1 demonstrate that the thermal design sensitivity 
calculated by the direct differentiation method is more accuracy than the results 
calculated by adjoint variable method. Moreover, by using the direct differentiation 
method, one can also get the time histories of the design derivatives of state variables 
as shown in Figs. 2 and 3. 

EXAMPLE 2: p-VERSION, TIME-DISCONTINIOUS FINITE ELEMENT 
APPROXIMATION 


This example will examine a more general heat conduction equation than the one 
presented in Eq. 1. The governing differential equation is given as 


du 3 3 3 
pc~r~ 2 S — 
dt ■=> j= 1 3fcc, 

u = f(x,t) 


k — 
* ax, 


= Q(x,t), in Qx(0,7°] 
on da u x(0,T] 


3 3 r l 

? X L, —[n,] = q,(x,t), 

'=•/=• dXj 
3 3 0M r i 

X2L,— [n^-hfu-T^), 
1=1 J= l ox, 


on dQ. q x(0,T] 
on 3Q A x(0,r] 


(36) 

(37) 

(38) 

(39) 


and 


u = g(x ) , 


in Q, at t = 0 


(40) 


where the temperature, «(x,t), is the only unknown. Eq. 36 represents an initial- 
boundary value problem with Eqs.37-39 being the temperature, the heat flux 



boundary condition, and the thermal convective condition; respectively, and Eq. 40, 
the initial conditions. It is assumed that the heat source, Q, the prescribed 
temperature, /j the flux, q % the thermal film coefficient, h, and the initial value, g, are 
all with proper regularity. In case of material nonlinearity, the specific heat, pc(u), the 
film coefficient, h(u), and the thermal conductivity, £y(u) are assumed to be functions 
of temperature. 


Let /= 0 on Qu- Otherwise, u in Eqs.36-40 can be replaced by u -f to achieve a 
homogenous boundary condition on As a result, the weak form of Eqs. 36-40 can 
be derived, based upon the Galerkin’s Method for an arbitrary function, w(x), as 


\(pc*^-w+ £ £ ky — ^ )dv+ \huwds 
l dt --U=1 dx t dXj £ 

= jQwdv+ J q s wds+ jhT^wds 
a dQ aa. 


(41) 


In the approach of time-discontinuous Galerkin’s method that is under 
development in this study, the time coordinate is treated as the same as the spatial 
coordinate. The time space is also discretized into elements or intervals. Focusing on 
one time interval, /„ = [t„.j, t „ ), the weak form, Eq. 41, can be extended to the entire 
product domain /„xi2 



n-i a 


du 33 du dw 

— -w+m. — — 

dt 1=1 /=• dx t dXj 


)dv ]dt + J7 jhuwds ]dt 

n - 1 dQ. h 


= j[ jQwdv + jq s wds+ jhT„wds]dt (42) 

n-1 £2 dQ. q 3 Q h 

where, w(x,t) is the testing function. Furthermore, to enforce the continuous 
requirement at the interfaces of time intervals, a weighted integral form of constraint, 
is appended to Eq. 42 

J (/x + u + - pc~u~)w + dv = 0 (43) 

a 


In this example, the temperature, u( x, t), is interpolated by a p- version hierarchical 
basis functions as 


u(x,t) =/ ( x,y,z,t )a 

=(**, ,y)®y(z)®m T a 


( 44 ) 


The symbol, ® , represents the outer tensor product operator and the vectors, yrznd 
0, represent collections of basis functions. Particularly, the through-thickness basis 
functions, yr, are made of the Legendre polynomials of the first kind 15 , the temporal 
basis functions, 0, the integrations of the same Legendre polynomials and the in-plane 
basis functions, chosen for triangular elements, as described by reference 21. 

It is assumed that in this study, the film coefficient h is set to zero and the 
relations between the other material properties, pc and Ky , and the temperature are 

defined in a tabulated form, presented as a result of experiments. Therefore, the 
material properties cannot be explicitly specified as functions of position and time as 
required by integration. An approximation is thus introduced to overcome such a 
difficulty. The standard Lagrange polynomials are used here for this purpose. 

The values of the material properties at the Lagrange points are taken from a 
given material table based upon the values of the temperature found at those points. 
The values of the material properties at elsewhere in an element are then obtained 
through interpolation. In this way, the material properties can be explicitly 
approximated as functions of position and time throughout the problem do main . 

As an example, the material property, say pc, can be interpolated in an element 
(/„ x^,) as 

pc(x,t) = (N c (x, y, z,t)) T £ 

= (N q ,( X ,y)®N cz (.z)®N ct (t)) T £ (45) 

= Z Z Z ( *0, (*, . ( Zj )®N ak (t i )) T S ciJk 

i j k 

where, 5 cijk is a component of the vector, £, which takes the value of the material 
properties found in the material table. 

If the material properties are interpolated linearly through the tabulated data, the 
value of the materials can be represented as, using pc as an example, 

^cijk = a yk + Pijk U ijk 

where a ijk and fi ljk are constants taken from the material table based upon the 
temperature value, u ijk , measured at the location (x t ,y t ,Zj ) and at time t k as 


“yk =zlka 

= (# , y, ) ® WiZj ) ® 6{t k )) T a 



SYSTEM EQUATION 


Once the interpolation functions are selected and the nonlinear material properties 
are approximated as explicit functions of position and time, one can proceed to 
integrate the terms in Eqs.42-43 to construct the equivalent matrices. The resultant 
finite element matrix equation for the time interval, /„ , can then be expressed as 

[c„ («„ )+ K n («») + k )k = 9n + W 'n-1 («„- 1 )]*»_! + (46) 


where the subscript, n, indicates that the associated matrix or vector is evaluated with 
functions defined in time interval, /„ . Note that each term in Eq.46 is corresponding 
to an integral in Eqs. 42-43, which can be evaluated based upon the interpolation 
functions of Eqs. 44-45. As an example, the capacitance matrix is given as 

C. - 

P 

The details of an elemental capacitance matrix, C a ^ is given as Eq. A.l in Appendix. 


The matrix equation derived above is based upon the discontinuous Galerkin’s 
method and will be solved in a time-marching fashion. In other words, the matrix 
equation of Eq. 46 will be solved for a time interval at a time, with a n as unknown 
and a„-j as known quantities. Because its nonlinear nature, Eq. 46 can be solved by 
the Newton-Raphson’s method, which leads to a recursive formula 


[c fl (ar , )+^fc' , )+<k' , )+^k' 1 )+jj < 1 )hi =-k 1 (47) 


and the solution is updated by 



i-l 

n 


+ A< 


(48) 


In the above equation, A a‘ n is the improvement of the solution and R‘ n 1 is the 
residual of the nonlinear equation at the i-l iteration, which is defined as 


R? n-rK 


Moreover, the derivative matrices, J c , Jt and J m are obtained by differentiating 
the coefficient matrices C„ , K„ and M* with respect to the unknown vector a„ , 
respectively. This is usually accomplished at the element level. As an example, the 
derivative matrix, J c , is detailed in Appendix as Eq. A.2. 


SENSITIVITY ANALYSIS 


Since the matrix equation, Eq. 46 , is in the form of a static problem, Ax=b. 
Differentiation of Eq. 46 with respect to a design variable, b, gives a sensitivity 
equation for a nonlinear transient heat transfer problem as, 

(C + K + M + +J C + / t + y„ + ).^- 


fdC dK dM + '] dM~ 
[db db db ) n db 


+ M~ 


fon-i , dq 
db db 


(49) 


The derivative matrices appearing in the right-hand side of Eq. 49 are new for 
sensitivity analysis. Generation of such matrices can be tedious and prone to 
mistakes. A typical derivative matrix, dCldb, is given in Appendix. Fortunately, 
because of the nature of approximation, those derivative matrices can be obtained in 
the same way as the original matrices themselves. This becomes evident by 
comparing Eqs. A1 and A3 in Appendix. Furthermore, it is noted that Eq. 49 is a 
da 

linear equation of — — , which enjoys the same left-hand side coefficient matrices as 
db 

Eq. 46 . Thus, the factored matrices saved from the converged analysis, can be reused 
here to solve the sensitivity equation. That makes the sensitivity analysis 
computationally efficient. 

NUMERICAL RESULTS 

An academic example is used here to verify the derived equations. The space 
domain of the 2-D problem is a 1x1 square, which is discretized into two triangular p- 
version finite elements, as shown in Fig. 4. The heat source term, Q, is specifically 
selected so that the solution of the heat transfer problem, Eq.36, is 

u t ( x , y,t ) = xy(x - l)(y - 1 )t 5 


which gives homogeneous boundary conditions and zero initial condition. Both the 
material properties, pc and isotropic k , are assigned by a linear temperature relation. 


pc,k = 20. + m 



The time step is set to be 1 second and the total operational time, t, is set to be 4 
seconds. 

In the numerical exercise, the orders of in-plane and the temporal polyno mials are 
selected to be 4 and 6, respectively, for temperature interpolation, Eq. 44, whereas the 
orders of in-plane and the temporal polynomials are 3 and 5, respectively, for material 
property interpolation, Eq. 45. The total Lagrange points for the in-plane material 
interpolation is 10 as marked in Fig. 4. Moreover, the error is measured by the L 2 - 
norm 15 as 



dxdy \dt 


1/2 


where e = u-u t and u is the numerical solution. 

At the end of the first time interval, r, = 1 second, the error is in the order of 10" 4 . 
The error is growing with the time. At t 2 = 2 second, the error grows to the order of 
10 1 . At f 4 =4, the error becomes 13.31. The major source of the error is expected 
from deficiency in the in-plane interpolation of materials. The material property is 
linear in temperature. Thus, the exact order of in-plane material interpolation has to 
be 4 which is higher than 3, the order used in the current study. Such error will be 
accumulated from one time interval to the next. The results of temperature 
distribution are shown in Fig. 5. 

Case 1 


The slope of the thermal conductivity-temperature relation is considered as a 
design variable. Since only the matrix K depends upon this special design variable, 
the right-hand side of Eq.49 can be greatly reduced to a single term - (dK / db) a. The 
results of the thermal sensitivity coefficients at times from 1 second to 4 seconds are 
shown in Fig. 6. Note that the thermal sensitivity coefficients is interpolated in the 
same way as the temperature; i.e.. 


du t . , da 

— = Y (x,t ) — 
db A db 


where da / db is obtained from Eq. 49. Comparing with the finite differencing, the 
errors of the thermal sensitivity coefficients are less than 1CT 4 for all the time 
intervals. Note that in this example, thermal sensitivity analysis takes only 18% of the 
time required for one thermal analysis. 


Case 2 


Since in this study, the material property is a distributed function, one may then 
assume that the square slab is made of various materials. In this particular case, the 
value of k at each Lagrange point is determined by it’s own material table. If the slope 
of each material table is considered as a design variable, there are 10 independent 
design variables in total. They are marked in Fig. 4. 

Figures 7 and 8 show the thermal sensitivity coefficients of du/dbj and du /dbg, 
where bj and bg are the slopes of the thermal conductivity-temperature relation at 
Lagrange points 1 and 9, respectively. The figures reveal that the design variable, bj, 
effects the change of temperature along the diagonal line, whereas the design 
variable, b^ does the area off the diagonal line more. Finally, all 10 thermal 
sensitivity coefficients of the temperature at the center point are collected and plotted 
out in Fig. 10. The picture indicates the degree of influence of individual design 
variable on the temperature at the center location at different time. 

CONCLUSIONS 

The paper uses two examples to demonstrate the derivation procedure for thermal 
sensitivity analysis. The continuous approach is used in the first example and the 
discrete approach is used in the second example. It is shown in the first example that 
the direct differentiation method can achieve better accuracy in thermal sensitivity 
analysis than the adjoint variable method. Several authors have similar observation. 
Furthermore, Example 1 shows that the direct results of the direct differentiation 
method, thermal derivatives of the temperature, are very useful in design. This 
particular advantage of the direct differentiation method is also demonstrated in 
Example 2. 

The second example presented here only represented an initial attempt to find the 
thermal sensitivity based upon the p-version time-discontinuous Galerkin’s method. 
Though construction of the matrix equation for thermal analysis is complicated, 
construction of that for thermal sensitivity analysis is rather simple. The resultant 
sensitivity equation is also demonstrated to be computationally efficient. However, 
more works are needed to develop error-control capability for thermal analysis and 
sensitivity analysis to ensure the quality of the p-version time-discontinuous 
Galerkin’s method. 
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Figure 1: Thermal Design Derivatives for Press Molding with Respect to the 

Mold Temperature. 
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Figure 2: Profiles of Thermal Design Derivatives of Temperature with Respect 

to the Mold Temperature. 



Figure 3: Profiles of Thermal Design Derivatives of the State of Cure with 
Respect to the Mold Temperature. 



Figure 4: Meshes and Lagrange Points for Interpolation of Material Property 
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Figure 6: Coefficients of Thermal Sensitivity with Respect to Homogeneous 
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Figure 7: Coefficients of Thermal Design Sensitivity with Respect to Thermal 

Conductivity at Material Lagrange Point 1; _ 
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Figure 8: Coefficients of Thermal Design Sensitivity with Respect to Thermal 

Conductivity at Material Lagrange Point 9; — . 
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Figure 9: Coefficients of Thermal Design Sensitivity with Respect to Thermal 
Conductivity; ; i =1 to 10. 
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Attachment II 

Fortran Source Code and Sample Input Data Files 

The files saved on the attached CD are listed below. 

1 . README.txt - Short Explanation of the attached FORTRAN Files 

2. schtmain6.f - the main source code based upon the p- version Galerkin’s 
method for a 3D nonlinear, transient heat conduction problem. 

3. input3dn6.dat - connectivity information for elements in time and 
through-thickness, the tabulated data for material properties 

4. joint_input.dat - connectivity information for elements in xy plane 

5. q3dt5zlxy4_pt6_pxy4_pz6nl.dat - input of thermal load vector which is 
obtained based upon the exact solution given by Eq.(20) 


